Skip to content

Conversation

@trontrytel
Copy link
Member

This PR fixes the latent heat bug discovered by @jbphyswx in our buoyancy gradients.

It would be great to rewrite parts of that file to get rid of the getter functions at the top and the bg_model, and to reduce the number of functions we define. But maybe we can do it in the next PRs.

get_qi_sat(thermo_params, bg_model) *
TD.latent_heat_sublim(thermo_params, t_sat)
) /
(get_ql_sat(thermo_params, bg_model) + get_qi_sat(thermo_params, bg_model))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add an eps to make sure this is not divided by zero?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good call. I thought I would be able to avoid it since we are only in that function if we have some non-zero cloud condensate. But I guess the numerics were not playing well with me - I actually need and sqrt(eps(FT)) in the denominator to make the implicit Rico single column test case work

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorry, I meant max(..., eps(FT)), but I can change it later. (or does max not work?)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did it change the results of any of the CI? I'm not super familiar with all the plots.
Also was wondering if it should converge to a real L_v at q=0 instead of 0, not that it should make a huge difference in the presence of any real cloud. I guess in equilibrium, condensate converges to 0 at the same time as the supersaturation, but that's less true in noneq.

sqrt(eps(Float64)) is only 1.49e-8 which is not unheard of for thin ice clouds, and if we ever ran Float32 simulations is 3e-4 which isn't insignificant (though I don't know if that's intended to be possible)

When I trialed fixing bg I just set to lh to Lv at ql+qi < eps(FT) -- I guess you could also just add +eps(FT)L_v or something similar to the numerator and denominator if you wanted it to be more numerically stable and continuous.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

prognostic edmf dycoms in ci changed. There are probably some other changes too. And we always run Float32. I agree we should do this more carefully. Could you work together with @trontrytel to make this better?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I have an idea. We can use TD.liquid_fraction here? That way Lh is always between Lv and Ls.

Copy link
Member

@szy21 szy21 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, just one comment.

@trontrytel trontrytel force-pushed the aj/buoyancy_gradients_lh branch from 7858c8f to 1a40a8c Compare November 20, 2025 23:44
@trontrytel trontrytel enabled auto-merge November 20, 2025 23:46
@trontrytel trontrytel added this pull request to the merge queue Nov 21, 2025
Merged via the queue into main with commit cec0868 Nov 21, 2025
19 checks passed
@trontrytel trontrytel deleted the aj/buoyancy_gradients_lh branch November 21, 2025 03:33
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants